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Abstract. A major advance in density-matrix renormalization group (DMRG) calculations has 
been achieved by the invention of highly efficient DMRG techniques for the simulation of real- 
time dynamics of strongly correlated quantum systems in one dimension. Starting from established 
linear-response techniques in DMRG and early attempts at real-time dynamics, we go on to review 
two current methods which both implement the idea of adapting the effective Hilbert space of 
DMRG to the quantum state evolving in time. We also give an outlook on extensions to finite 
temperature calculations. 



INTRODUCTION 

The physics of strongly correlated quantum systems continues to pose major challenges 
in experimental and theoretical physics. While both experiment and theory have focused 
on static, thermodynamic or at most linear-response quantities in the past, recently ques- 
tions which explicitly involve the out-of-equilibrium time-dependence of such quantum 
systems have come to the foreground. These questions arise in the context of transport 
far from equilibrium or of decoherence, particularly as the size of devices continues to 
shrink towards the atomic scale. However, perhaps the most striking example for this 
is provided by the progress in preparing dilute ultracold bosonic and also fermionic al- 
kali gases. Subjected to an optical lattice, these gases are arguably the purest realization 
of the typical model Hamiltonians of strong correlation physics, such as the Hubbard 
model[ 1, 2]. More importantly, the interaction parameters can be tuned experimentally 
on quantum mechanically relevant time-scales over a huge range, while being known 
precisely from microscopic calculations. From a theoretician's point of view, this sit- 
uation is almost ideal, and has stimulated great interest in the development of time- 
dependent methods. 

For linear response, exact diagonalization can provide detailed results for small sys- 
tems, and quantum Monte Carlo can provide coarse resolution results after analytic con- 
tinuation from imaginary time, for systems without the sign problem. Outside the linear 
response regime, almost the only tool available has been the diagonalization of very 
small clusters. 

In this review, the emphasis is on recent extensions of the density-matrix renormal- 
ization group method (DMRG) into the real-time domain which make it the 



currently most powerful method for such problems. Following up on early attempts to 
extend DMRG to real-time, input from quantum information theory has led to the for- 
mulation of two DMRG algorithms for real-time evolutions. We set out with a reminder 
about the basic ideas of DMRG, review linear response calculations with DMRG and 
move on to early attempts in the time domain. We then explain how the TEBD algo- 
rithm of Vidal [6] beautifully reflects fundamental structures of DMRG and hence can 
be easily used to extend DMRG to the time-domain |7, 8]. However, this approach has 
shortcomings for longer-ranged interactions, which can be circumvented in yet another 
modification of DMRG at some cost in efficiency The range and power of both 
methods are discussed based on "real-life" applications. 



BASIC IDEAS OF DMRG 

Several good descriptions of the DMRG algorithms exist in the literature . Rather 

than repeat these descriptions, here we summarize the most important ideas of DMRG. 

The first key idea is the description of a collection of sites, or block, in terms of a 
limited set of basis states and operator matrices between those states. These states and 
matrices are defined by a set of basis transformations as sites are successively added 
to the block. This representation is due to Wilson and is a key feature of his numerical 
renormalization group (NRG^ lllOll . Let the block at the beginning of step i be described 
by a set of states In this step site i (states {|c7)}) is added to the block. The new 
states describing the larger block { | } are given by 

|0=lAf,,[a]|/)®|a). 

a,i 

The number of states is kept approximately constant at m, so the set of states 
is incomplete. If the states {|/')} were described in detail in terms of the sites, the 
computational effort would grow exponentially despite the truncation to a constant 
number of states. Instead, the m x m matrices for the Hamiltonian and other operators 
give all the detail needed to construct the Hamiltonian at larger length scales. These 
matrices are transformed to the new basis at each step using the transformation matrices 
A. In this way, the computational effort remains constant as 0{m^). 

The second key idea of DMRG is to choose the states to keep as eigenstates of 
the reduced density matrix (RDM) of the block. In Wilson's NRG approach, one kept 
the lowest energy eigenstates of the block Hamiltonian. This choice works for the 
special Hamiltonians devised by Wilson for impurity systems, but fails for more general 
lattice systems. Using the density matrix eigenstates can be shown to be optimal for 
reproducing the wavefunction as well as the RDM. Let the block have states \i) and the 
rest of the system, referred to as the environment, have states Then the wavefunction 
of the whole system is written as 

\¥)=Y.¥ij\i)^\j) 

ij 
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FIGURE 1. The standard DMRG superblock configuration, in which the left central site is being added 
onto the left block. 



and the coefficients of tlie RDM p are 

j 

The eigenvalues of p are the probabilities of the block being in the corresponding 
eigenstate, and if the probability is neglible, the eigenstate can be omitted from the basis. 
The RDM can also be built by summing over several states 1 with arbitrary weights, 
representing the probability of each in a mixed state of the system. In this case each 
I Xf/^) is said to be targetted, an important concept for time-dependent DMRG. 

The RDM depends on the enviroment through 1 1//) . In DMRG, both the block and 
the environment are described approximately using basis sets of size ~ m. This leads 
to the third key idea of DMRG, the idea of sweeping back and forth to produce a 
self-consistent, accurate representation for both parts of the system. In Fig. 1 we show 
the most common superblock configuration in DMRG. In the sweeping procedure the 
dividing line between the left and right block moves back and forth between the ends of 
the system. The block which is growing is treated as the system block; the other is the 
enviroment, with the roles reversed when the direction is reversed. The system block in 
each case obtains an improved basis during the sweep. In Wilson's NRG, there was no 
sweeping and no feedback from the low energy scales to the high energy scales. For a 
general lattice system, not divided by energy scales, feedback is necessary. 

If we trace back through the basis set transformations which led to the basis of a block, 
one obtains an explicit representation of the states 

\h)= E [A'[ai]A^[a2]...A''[ae]Uai...ae) 

ai...ae 

where A Ms a vector for each Oi and the rest of the A's are matrices. We can write the 
basis states for the right block similarly. Alternatively, we can let the dividing line be all 
the way on the right end of the system, in which case we can write the following matrix 
product expression for | i/a) : 

\W)= L A\ai]A^[o2]...A'^[aL]\a,...aL) 

where A^[aL] is a vector for each value of Oi, so that the product of the A's is a scalar. 
The matrix product representation for 1 1//) was first developed in a DMRG context by 
Ostlund and Rommer lUlll . but its usefulness was not widely appreciated at the time. 
More recently, the matrix product representation has become very important as a route 
to improve the capabilities of and to generalize DMRG [.12.1 . 
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FIGURE 2. Diagrams for the matrix product representation of the wavefunction and an operator on a 
block. 



A useful diagrammatic form of the matrix product representation (v^ is illustrated in 
Fig. 2. The upper diagram represents the wavefunction. Here each intersection of lines 
is associated with a site and represents a matrix (here the A^[a^]), or more generally, a 
tensor. The vectors and are reprecented by right-angle segments. Interior line 
segments represent indices which are summed over, whereas the ends of segments 
sticking out of the figure represent external indices labeling states. As another example, 
let an operator O be defined as a product of site operators on the left block, Oi. . .Og, 
where many of the Oi may be identity operators. The lower diagram represents the 
matrix for this operator in the basis of the left block. Each of the two lines sticking 
out on the right represents indices running over the m states of the left block. 



DMRG DYNAMICS 

While the original DMRG algorithm seems to be limited to the calculation of equi- 
librium properties such as ground state correlations, it can in fact be extended to linear 
response quantities. For some operator A, we define a (time-dependent) Green's function 
at r = in the Heisenberg picture by 

iGA{t'-t) = {0\A\t')Ait)\0) (1) 

with t' >t for a time-independent Hamiltonian^. Going to frequency-space, the Green's 
function reads 

GA{(0 + iil) = {0\P-—-L_ ^A|0), (2) 

Eq + co + iti —H 

where rj is some positive number to be taken to zero at the end. We may also use the 
spectral or Lehmann representation of correlations in the eigenbasis of H, 

CAico) = Y,\{n\A\0)\^5{co + Eo-En). (3) 

n 

Ca{co) is related to Ga{co + irj) as 

Ca{o})= lim -—ImGAico + irj). (4) 
77^0+ 7t 



The role of rj in DMRG calculations is threefold: First, it ensures causality in Eq. 
©. Second, it introduces a finite lifetime t oc l/rj to excitations. Third, 77 provides 
a Lorentzian broadening of Ca((u), 

C.(a. + i^) = i/</«'C.(a,')j^^-^L_, (5) 

which serves either to broaden the numerically obtained discrete spectrum of finite 
systems into some "thermodynamic limit" behavior or to broaden analytical results for 
Ca for comparison to numerical spectra where 77 > 0. 

Most DMRG approaches to dynamical correlations center on the evaluation of Eq. 
The first, which we refer to as Lanczos vector dynamics, has been pioneered by 
Hallberg lUsll . and calculates highly time-efficient, but comparatively rough approxima- 
tions to dynamical quantities adopting the Balseiro-Gagliano method iTull to DMRG. 
The second class of approaches, including both the correction vector method I 15l IT6I1 
and DDMRG (dynamical DMRG) [ 17J, is also based on pre-DMRG techniques, but is 
both much more precise and numerically much more expensive. 

Boundary effects due to DMRG-typical open boundary conditions can be treated in 
various ways; two situations should be distinguished. If A acts locally, such as in the 
calculation of an optical conductivity, one may exploit that finite rj exponentially sup- 
presses excitations 1 17]. As they travel at some speed c through the system, a thermo- 
dynamic limit L — > 00 first, rj ^ second may be taken consistently as a single limit 
with 77 = c/L. For the calculation of dynamical structure functions such as obtained in 
elastic neutron scattering, A is a spatially delocalized Fourier transform, and another ap- 
proach must be taken. The open boundaries introduce both genuine edge effects and a 
hard cut to the wave functions of excited states in real space, leading to a large spread 
in momentum space. To limit bandwidth in momentum space, filtering by modifying 
A{x) A{x)f{x) is necessary. The filtering function f{x) should be narrow in momen- 
tum space and broad in real space, while simultaneously strictly excluding edge sites. 
For a detailed discussion of such filters, see [16]. 



Continued fraction dynamics 

The technique of continued fraction dynamics has first been exploited by Gagliano 
and Balseiro [14] in the framework of exact ground state diagonalization. Obviously, 
the calculation of Green's functions as in Eq. © involves the inversion of H (or more 
precisely, £0 + ft) + atypically very large sparse hermitian matrix. This inversion 

is carried out in two, at least formally, exact steps. First, an iterative basis transformation 
taking ^ to a tridiagonal form is carried out. Second, this tridiagonal matrix is then 
inverted, allowing the evaluation of Eq. @. 

Let us call the diagonal elements of H in the tridiagonal form a„ and the subdiagonal 
elements Z?^. The coefficients a„, are obtained as the Schmidt-Gram coeffcients in 
the generation of a Krylov subspace of unnormalized states starting from some arbitrary 
state, which we take to be the excited state A |0): 



\fn+l)=H\fn)-an\fn)-bl\fn-l). 



(6) 



with l/o) =A|0), and 



^ {fn\H\fn) 2 ^ {fn-i\H\fn) ^ {fn\fn) 

(fnlfn) ' " {fn-l\fn-l) (/n-ll/n-l)" ^ ^ 

The global orthogonality of the states (at least in formal mathematics) and the 
tridiagonality of the new representation (i.e. {fi\H\fj) = for \i — j\ > 1) follow by 
induction. It can then be shown quite easily by an expansion of determinants that the 
inversion of Eq + co + irj — H leads to a continued fraction such that the Green's function 
Ga reads 

<°l^*'l°> ■ (8) 

z-ao '-a- 

where z = EQ-\-(0 + ir\. This expression can now be evaluated numerically, giving access 
to dynamical correlations. Alternatively, one may also exploit that upon normalization of 
the Lanczos vectors and accompanying rescaling of the a„ and b^, the Hamiltonian is 
iteratively transformed into a tridiagonal form in a new approximate orthonormal basis. 
Transforming the basis { } by a diagonalization of the tridiagonal Hamiltonian matrix 
to the approximate energy eigenbasis of H, {\n)} with eigenenergies £„, the Green's 
function can be written within this approximation as 

GA((0 + i77)= (9) 

where the sum runs over all approximate eigenstates. The dynamical correlation function 
is then given by 

where the matrix elements in the numerator are simply the |/o) expansion coefficients 
of the approximate eigenstates \n). 

In practice, several limitations occur. The iterative generation of the coefficients a„, 
is equivalent to a Lanczos diagonalization of H with starting vector A |0). Typically, 
the convergence of the lowest eigenvalue of the transformed tridiagonal Hamiltonian to 
the ground state eigenvalue of H will have happened after n ~ 0(10^) iteration steps 
for standard model Hamiltonians. Lanczos convergence is however accompanied by 
numerical loss of global orthogonality which computationally is ensured only locally, 
invalidating the inversion procedure. With A|0) as starting vector, convergence will be 
fast if A 1 0) is a long-lived excitation (close to an eigenstate) such as would be the case if 
the excitation is part of an excitation band; this will typically not be the case if it is part 
of an excitation continuuum. Moreover, H itself is not exact, and its repeated application 
generates further errors. 

As an example for the excellent performance of this method, one may consider the 
isotropic spin-1 Heisenberg chain, where the single magnon line is shown in Fig. El Ex- 
act diagonalization, quantum Monte Carlo and DMRG are in excellent agreement, with 
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FIGURE 3. Single magnon line of the 5=1 Heisenberg AFM from exact diagonalization, quantum 
Monte Carlo and DMRG for various system sizes and boundary conditions. From LI 6,1 . 




FIGURE 4. Spectral weight S^{q — 7Z, Co) of the 5=1/2 Heisenberg AFM from Lanczos vector and 
correction vector DMRG. Nl indicates the number of target states; M = 256. Note that spectral weight 
times CO is shown. From L16.I . 



the exception of the region ^ ^ 0, where the single-magnon band has disappeared into 
a two magnon continuum. Here Lanczos vector dynamics does a poor job reproducing 
the peak of the spectral function just above the bottom of the continuum, which has a 
gap of twice the Haldane gap Ah at ^ = 0. 

The intuition that excitation continua are badly approximated by a sum over some 
0(10^) effective excited states is further corroborated by considering the spectral weight 
function S^{q= 7r,(0) [use A = 5+ in Eq. ©] for a spin- 1/2 Heisenberg antiferromagnet. 
As shown in Fig.lH Lanczos vector dynamics roughly catches the right spectral weight, 
including the 1 /ft) divergence, as can be seen from the essentially exact correction vector 
curve, but no convergent behavior can be observed upon an increase of the number 
of targeted vectors. The very fast Lanczos vector method is thus certainly useful to 
get a quick overview of spectra, but not suited to detailed quantitative calculations of 
excitation continua, only excitation bands. 



Correction vector dynamics 



Even before the advent of DMRG, another way to obtaining more precise spectral 
functions had been proposed in [18]; it was first applied using DMRG in ifTsil and lH^ . 
After preselection of a fixed frequency CO one may introduce a correction vector 

k(« + ir])) = TT— ^A|0), (11) 
Eq + co + iT} —H 

which, if known, allows for trivial calculation of the Green's function and hence the 
spectral function at this particular frequency: 

GA{co + i'n) = {A\c{co + i'n)). (12) 

The correction vector itself is obtained by solving the large sparse linear equation system 
given by 

{Eo + co + iri-H)\c{a + ir]))=A\0). (13) 

To actually solve this nonhermitean equation system, the current procedure is to split the 
correction vector into real and imaginary part, to solve the hermitean equation for the 
imaginary part and exploit the relationship to the real part: 

[{Eo + (o-H)^ + ri^]lm\c{Q} + iri)) = -riA\0) (14) 

Re\c{(0 + iri)) = ^ Im|c((0 + i7])) (15) 

The standard method to solve a large sparse linear equation system is the conjugate- 
gradient method, which effectively generates a Krylov space as does the Lanczos al- 
gorithm. The main effor in this method is to provide ^^Im|c). Two remarks are in 
order. The reduced basis representation of is obtained by squaring the effective 
Hamiltonian generated by DMRG. This approximation is found to work extremely well 
as long as both real and imaginary part of the correction vector are included as tar- 
get vectors: While the real part is not needed for the evaluation of spectral functions, 
{Eq + CO — ^)Im|c) ~ Re|c) due to Eq. (fT51) : and targeting Re|c) ensures minimal trun- 
cation errors in ^Im|c). The fundamental drawback of using a squared Hamiltonian is 
that for all iterative eigenvalue or equation solvers the speed of convergence is deter- 
mined by the matrix condition number which drastically deteriorates by the squaring of 
a matrix. This might be avoided by using biconjugate or conjugate symmetric equation 
solvers for Eq. ([T3t directly. 

Alternatively, there is a reformulation of the correction vector method in terms of 
a minimization principle, which has been called "dynamical DMRG" [17]. While the 
fundamental approach remains unchanged, the large sparse equation system is replaced 
by a minimization of the functional 

WA,^{co,xt/)= (16) 
{\ir\{Eo + co-H)^ + ri^\\ir) + n{A\\ir) + n{xif\A). 
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FIGURE 5. Spectral weight of the 5 = 1 /2 Heisenberg AFM from correction vector DMRG. M = 
256 states kept. Spectral weights have been calculated for co-intervals starting from various anchoring 
frequencies for the correction vector. From flfill . 

At the minimum, the minimizing state is 

IWin) =Im|c((D + ir])). (17) 

Even more importantly, the value of the functional itself is 

WA,Tji(o,\j/) = -7triCA{(0 + iri), (18) 

such that for the calculation of the spectral function it is not necessary to explicitly use 
the correction vector. In the simplest form of the correction vector method, the density 
matrix is formed from targeting four states, |0), A|0), Im|c(£0 + i77)) and Re|c(ft) + iT])). 

As has been shown in [16], it is not necessary to calculate a very dense set of 
correction vectors in co-space to obtain the spectral function for an entire frequency 
interval, assuming that the finite convergence factor rj ensures that an entire range of 
energies of width ^ 7] is described quite well by the correction vector. It was found 
that best results are obtained for a two-correction vector approach where two correction 
vectors are calculated and targeted for two frequencies cOi, (O2 = + Aft) and the 
spectral function is obtained for the interval [oji , coz] using the Lanczos method for 
the approximate Hamiltonian produced by this targeting scheme. This method is, for 
example, able to provide a high precision result for the spinon continuum in the 5=1/2 
Heisenberg chain where standard Lanczos dynamics fails (Fig. 

EARLY ATTEMPTS AT TIME-EVOLUTION 

Even though the methods described in the previous section provide high-quality linear- 
response quantities, they fail in truly out-of-equilibrium situations or for time-dependent 
Hamiltonians; where they work, they are very time-consuming. It has therefore been of 
high interest to find DMRG approaches dealing with state evolution in real-time. 

To see the advantages of such an approach, consider the following. Essentially all 
physical quantities of interest involving time can be reduced to the calculation of either 



equal-time n-point correlators such as the (1 -point) density 



= m)\niW{t)) = (v/k'^V-'^» (19) 
or unequal-time n-point correlators such as the (2-point) real-time Green's function 

G„(0 = {W]i.t)cj{Q)W) = {Mr\e+-^^'c]e--^"'cj\xir). (20) 

This expression can be cast in a form very close to Eq. ( fT^ by introducing |0) = Cj|i//) 
such that the desired correlator is then simply given as an equal-time matrix element 
between two time-evolved states, 

G,j{t) = {Mrit)\cJ\(l>it)). (21) 

If both |v^(?)) and \(j>(t)) can be calculated, a very appealing feature of this approach 
is that Gij{t) can be evaluated in a single calculation for all i and t as time proceeds. 
Frequency-momentum space is then reached by a double Fourier transformation. Obvi- 
ously, finite system-sizes and edge effects as well as algorithmic constraints will impose 
physical constraints on the largest times and distances \i — j\ or minimal frequency and 
wave vectors resolutions accessible. Nevertheless, this approach might emerge as a very 
attractive alternative to the current very time-consuming calculations of G{k, co) using 
the dynamical DMRG| 16, 17]. 

The fundamental difficulty of obtaining the above correlators becomes obvious if we 
examine the time-evolution of the quantum state 1 = 0) ) under the action of some 
(for simplicity) time-independent Hamiltonian i?| = En\^^n)■ If the eigenstates | ^/n) 
are known, expanding | = 0)) =Y.n Cn\Wn) leads to the well-known time evolution 

IV^(O) = £c„exp(-i£„f)| v/„), (22) 

n 

where the modulus of the expansion coefficients of | ^^{t)) is time-independent. A sensi- 
ble Hilbert space truncation is given by a projection onto the large-modulus eigenstates. 
In strongly correlated systems, however, we usually have no good knowledge of the 
eigenstates. Instead, one uses some orthonormal basis with unknown eigenbasis expan- 
sion, \k) = Y.n^kn\¥n)- The time evolution of the state = 0)) = Y.kdk{^)\k) then 
reads 

IV^(O) = E (Y.dk{0)akne-'^A \k) = £4(01-^), (23) 

k ^ n / k 

where the modulus of the expansion coefficients di^(t) is time-dependent. For a general 
orthonormal basis, Hilbert space truncation at one fixed time (i.e. t = 0) will therefore 
not ensure a reliable approximation of the time evolution. Also, energy dijferences 
matter in time evolution due to the phase factors e^i(^"^^n')'' in \dk{t)\^. Thus, a good 
approximation to the low-energy Hamiltonian alone (as provided by DMRG) is of 
limited use. 

Static time-dependent DMRG. Cazalilla and MarstonfT^ were the first to exploit 
DMRG to systematically calculate time-dependent quantum many-body effects. They 



studied a time-dependent Hamiltonian H{t) = H{0) + V{t), where V{t) encodes the 
time-dependent part of the Hamiltonian. After applying a standard DMRG calculation 
to the Hamiltonian H{t = 0), the time-dependent Schrodinger equation was numerically 
integrated forward in time. The effective Hamiltonian in the reduced Hilbert space 
was built as ^eff(0 = ^eff(O) +K;ff(0' where ^eff(O) was taken as the last superblock 
Hamiltonian approximating H{0). VefflO as an approximation to V was built using 
the representations of operators in the final block bases. The initial condition was 
obviously to take |V^(0)) as the ground state obtained by the preliminary DMRG run. 
This procedure amounts to working within a static reduced Hilbert space, namely that 
optimal at ? = 0, and projecting all wave functions and operators onto it. 

In this approach the hope is that an effective Hamiltonian obtained by targeting the 
ground state of the t = Hamiltonian is capable to catch the states that will be visited 
by the time-dependent Hamiltonian during time evolution. This approach must however 
break down after relatively short times as the full Hilbert space is explored, as became 
quickly obvious. 

Dynamic time-dependent DMRG. Several attempts have been made to improve on 
static time-dependent DMRG by enlarging the reduced Hilbert space using information 
on the time-evolution, such that the time-evolving state has large support on that dynamic 
Hilbert space for longer times. Whatever procedure for enlargement is used, the problem 
remains that the number of DMRG states m grows with the desired simulation time as 
they have to encode more and more different physical states. As calculation time scales 
as m^, this type of approach will meet its limitations somewhat later in time. 

All enlargement procedures rest on the ability of DMRG to describe - at some 
numerical expense - small sets of states ("target states") very well instead of just one. 

The simplest approach is to target the set {| V^,)} = {\^f{Q)),H\\if{Q)),H^\\if{0)),...}. 
Alternatively, one might consider the Krylov vectors formed from this set. Results 
improve, but not decisively. 

A much more time-consuming, but also much better performing approach has been 
demonstrated by Luo, Xiang and Wang lEol] . They use a density matrix that is given by a 
superposition of states | V^(?, ) ) at various times of the evolution, p = E/^o I ) ( ^i^i) I 
with E Of; = 1 for the determination of the reduced Hilbert space. Of course, these states 
are not known initially; it was proposed by them to start within the framework of infinite- 
system DMRG from a small DMRG system and evolve it in time. For a very small 
system this procedure is exact. For this system size, the state vectors \ ^f{ti)) are used to 
form the density matrix. This density matrix then determines the reduced Hilbert space 
for the next larger system, taking into account how time-evolution explores the Hilbert 
space for the smaller system. One then moves on to the next larger DMRG system where 
the procedure is repeated. This is of course very time-consuming. 

Schmitteckert[21] has computed the transport through a small interacting nanostruc- 
ture using an Hilbert space enlarging approach, based on the time evolution operator. To 
this end, he splits the problem into two parts: By obtaining a relatively large number of 
low-lying eigenstates exactly (within time-independent DMRG precision), one can cal- 
culate their time evolution exactly. For the subspace orthogonal to these eigenstates, he 
implements the matrix exponential | + A?)) = exp(— ii/A?)|i//(?)) using the Krylov 
subspace approximation. For any block-site configuration during sweeping, he evolves 



the state in time, obtaining \yf{ti)) at fixed times tj. These are targeted in the density 
matrix, such that upon sweeping forth and back a Hilbert space suitable to describe all 
of them at good precision should be obtained. For numerical efficiency, he carries out 
this procedure to convergence for some small time, which is then increased upon sweep- 
ing, bringing more and more states \ ^l^{ti)) into the density matrix. Again, this is a very 
time-consuming approach. 



TIME-EVOLVING BLOCK DECIMATION 

Decisive progress came from an unexpected corner, namely quantum information theory, 
when Vidal proposed an algorithm for simulating quantum time evolutions of one- 
dimensional systems efficiently on a classical computer flllS]- His algorithm, known as 
TEBD [time-evolving block decimation] algorithm, is based on matrix product states iEbI 
I24] : as it turned out, it is so closely linked to DMRG concepts, that his ideas could be 
implemented easily into DMRG, leading to an adaptive time-dependent DMRG, where 
the DMRG state space adapts itself in time to the time-evolving quantum state. In this 
section, we will explain his algorithm. 

A useful concept is that of a Schmidt decomposition: Consider a quantum state 
1 1^) = Y,ij Wij 1 ® I j) introduced before, with A^'^ states | /) and A^"^ states | j) . Assuming 
without loss of generality A^'^ > A^^, we form the {N^ x A^^) -dimensional matrix A with 
Aij = \ifij. Singular value decomposition guarantees A = UDV^ , where U is (A^"^ x A^^)- 
dimensional with orthonormal columns, D is a (A^^ x A^^) -dimensional diagonal matrix 
with non-negative entries Daa = ^/w^, and is a (A^^^ x A^^) -dimensional unitary 
matrix; 1 1//) can be written as 

pjS 

IV^) ^ mt^^«v^Vi-|0|j) (24) 

;=1 a=l j=l 

NE / j^S \ / ^E \ 

The orthonormality properties of U and ensure that \w^^) = Y^iUia\i) and \w^) — 
Y^jVja\j) form orthonormal bases of system and environment respectively, in which the 
Schmidt decomposition 

^'Schmidt 

|V^)= £ V^ka)ka) (25) 

holds. N^N^ coefficients are reduced to A/^schmidt < non-zero coefficients ^/w^, 
wi > W2 > > ■ ■ ■• Relaxing the assumption A'^'^ > A^^, one has 



A^Schmidt<min(A^^,A^^). 



(26) 
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FIGURE 6. Bipartitioning by cutting bond I between sites I and Z + 1 . 



Upon tracing out environment or system the reduced density matrices for system and 
environment are found to be 

^Schmidt ^Schmidt 

P5= £ Wa\w%){w%\; pE= ^a\wl){wl\. (27) 

a a 

DMRG reduced density matrix analysis and the Schmidt decomposition therefore yield 
exactly the same information. This fact was understood from the very beginning of 
DMRG, although we had not heard the term "Schmidt decomposition". In fact, the 
singular value decomposition representation of the wavefunction was understood before 
the density matrix representation. 

Let us now formulate the TEBD simulation algorithm. In the original exposition 
of the algorithm [221. one starts from a representation of a quantum state — 

LfTi...fTi V^c7i ctl|o'i---Ol) where the coefficients for the states are decomposed as a 

product of tensors, 

V^cT, ,...,a, = E T'a, [cTi] Tl^ [CT2] A^^F^ [a,] ■ ■ ■ r^^_, [at] . (28) 

Oil, 

It is of no immediate concern to us how the T and A tensors are constructed explicitly 
for a given physical situation. Let us assume that they have been determined such that 
they approximate the true wave function close to the optimum obtainable within the class 
of wave functions having such coefficients; this is indeed possible as will be discussed 
below. There are, in fact, two ways of doing it: within the framework of DMRG, or by 
a continuous imaginary time evolution from some simple product state, as discussed in 
Ref. 10]. 

The ansatz can be visualized: the (diagonal) tensors X\i= 1, . . . ,L — 1 are associated 
with the bonds i, whereas F', / = 2, . . . , L — 1 links (transfers) from bond i to bond / — 1 
across site i. Note that at the boundaries (z = 1 , L) the structure of the T is modified. The 
sums run over m states | a,) living in auxiliary state spaces on bond i. A priori, these 
states have no physical meaning. 

The r and A tensors are constructed such that for an arbitrary cut of the system into a 
part Si of length I and a part / of length L — ldX bond /, the Schmidt decomposition 
for this bipartite splitting reads 

IV^)=I<l4)l4r), (29) 

a, 

with 

W%)= E I rJ,JcTi]Ai---r'„^_^«Ja,]|ai)®---®|cTz), (30) 

..,«/_! CTi,...,CT; 



and 

|a/+i)®---®|aL), (31) 

where is normalized and the sets of {Iwa,)} and orthonormal. This 

implies, for example, that 

I«)' = l- (32) 

We can see that (leaving aside normalization considerations for the moment, see 
10]) this representation may be expressed as a matrix product state if we choose for 
A'[c7,-]=L«,^A'^^[a,-]|a)(/3| 

<p[<yi]=Kp[^-Wp^ (33) 

except for z = 1 , and i = L, where expressions are slightly modified. 

Let us now consider the time evolution for a typical (possibly time-dependent) Hamil- 
tonian with nearest-neighbor interactions: 

^ = EA'+i+ E Gj,j+u (34) 

i odd j even 

,+1 and Gjj+i the local Hamiltonians on the odd bonds linking i and z -|- 1, and the 
even bonds linking j and j +1. While all F and G terms commute among each other, F 
and G terms do in general not commute if they share one site. Then the time evolution 
operator may be approximately represented by a (first order) Trotter expansion as 

^-iH8t _ Yl e-iA/+i5? Yl e-'^J J+i^' + ff{5t^), (35) 

i odd j even 

and the time evolution of the state can be computed by repeated application of the 
two-site time evolution operators exp{—iGjj+i5t) and exp(— i/) ,_|_i5?). This is a well- 
known procedure in particular in Quantum Monte Carlo[25] where it serves to carry out 
imaginary time evolutions (e.g. checkerboard decomposition). 
The TEBD simulation algorithm now runs as follow sll5 l22l] : 

1. Perform the following two steps for all even bonds (the order does not matter): 

(i) Apply exp(— iG/,/+i5?) to |v^(0)- For each local time update, a new wave 
function is obtained. The number of degrees of freedom on the "active" bond 
thereby increases, as will be detailed below. 

(ii) Carry out a Schmidt decomposition cutting this bond and retain as in DMRG 
only those m degrees of freedom with the highest weight in the decomposition. 

2. Repeat this two-step procedure for all odd bonds, applying exp(— iF/ /_|_i5?). 

3. This completes one Trotter time step. One may now evaluate expectation values at 
selected time steps, and continues the algorithm from step 1. 



Let us now consider some computational details. 

(i) Consider a local time evolution operator acting on bond /, i.e. sites / and / + 1, for a 
state 1 1//) . The Schmidt decomposition of 1 1//) after partitioning by cutting bond / reads 

M 

\W)= K|wg)|w^r')- (36) 

0!/ = l 

Using Eqs. (l30b and (ISTT) . we find after expanding into and \oi), and 

similarly for Iw^J""'), 

(37) 

We note, that this has the form of a typical DMRG state for two blocks and two sites 

IV^) = L E E L V^™/-ICT,CT,+Im,+JW^,_,)|CT/)|C7/+I)|w^,^j). (38) 

The local time evolution operator on site 1,1+ I can be expanded as 

Ui,i^i= I I u''J;Jlll\alaU,){aiai+,\ (39) 

and generates | \]/') = | \]/)- 
where 

(ii) Now a new Schmidt decomposition identical to that in DMRG can be carried out for 

cutting once again bond Z, there are now mngite states in each part of the system, 
leading to 

w)= I <i4)i^«r')- (42) 

0!/ = l 

In general the states and coefficients of the decomposition will have changed compared 
to the decomposition (l36t previous to the time evolution, and hence they are adaptive. 
We indicate this by introducing a tilde for these states and coefficients. As in DMRG, 
if there are more than m non-zero eigenvalues, we now choose the m eigenvectors 
corresponding to the largest XL to use in these expressions. The error in the final state 
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FIGURE 7. Finite-system DMRG algorithm. Block growth and shrinkage. For the adaptive time- 
dependent DMRG, replace ground state optimization by local time evolution. 



produced as a result is proportional to the sum of the magnitudes of the discarded 
eigenvalues. After normalization, to allow for the discarded weight, the state reads 

m 

\V)= (43) 

Note again that the states and coefficients in this superposition are in general different 
from those in Eq. (l36b : we have now dropped the tildes again, as this superposition will 
be the starting point for the next time evolution (state adaption) step. 

The key point about the TEBD simulation algorithm is that a DMRG-style truncation 
to keep the most relevant density matrix eigenstates (or the maximum amount of entan- 
glement) is carried out at each time step. This is in contrast to previous time-dependent 
DMRG methods, where the basis states were chosen before the time evolution, and did 
not "adapt" to optimally represent the state at each instant of time. 



ADAPTIVE TIME-DEPENDENT DMRG 

DMRG generates position-dependent m x m matrix-product states as block states for a 
reduced Hilbert space of m states; the auxiliary state space to a bond is given by the 
Hilbert space of the block at whose end the bond sits. This physical meaning attached 
to the auxiliary state spaces implies that they carry good quantum numbers for all block 
sizes. The big advantage is that using good quantum numbers allows us to exclude a large 
amount of wave function coefficients as being 0, drastically speeding up all calculations 
by at least one, and often two orders of magnitude. 

The effect of the finite-system DMRG algorithm[4] is now to shift the two free sites 
through the chain, growing and shrinking the blocks S and E as illustrated in Fig. |7] At 
each step, the ground state is redetermined and a new Schmidt decomposition carried 
out in which the system is cut between the two free sites, leading to a new truncation 
and new reduced basis transformations (2 matrices A adjacent to this bond). 

As the actual decomposition and truncation procedure in DMRG and the TEBD 
simulation algorithm are identical, one can use the finite-system algorithm to carry 
out the sequence of local time evolutions (instead of, or after, optimizing the ground 



state), thus constructing by Schmidt decomposition and truncation new block states best 
adapted to a state at any given point in the time evolution (hence adaptive block states) as 
in the TEBD algorithm, while maintaining the computational efficiency of DMRG. To 
do this, one needs not only all reduced basis transformations, but also the wave function 
I in a two-block two-site configuration such that the bond that is currently updated 
consists of the two free sites. This implies that has to be transformed between 
different configurations. In finite-system DMRG such a transformation, which was first 
implemented by White lE9l ("state prediction") is routinely used to predict the outcome 
of large sparse matrix diagonalizations, which no longer occur during time evolution. 
Here, it merely serves as a basis transformation. 

The adaptive time-dependent DMRG algorithm which incorporates the TEBD simu- 
lation algorithm in the DMRG framework is therefore set up as follows: 

0. Set up a conventional finite-system DMRG algorithm with state prediction using 
the Hamiltonian at time t = 0, H{0), to determine the ground state of some system 
of length L using effective block Hilbert spaces of dimension M. At the end of this 
stage of the algorithm, we have for blocks of all sizes / reduced orthonormal bases 
spanned by states |m/), which are characterized by good quantum numbers. Also, 
we have all reduced basis transformations, corresponding to the matrices A. 

1 . For each Trotter time step, use the finite-system DMRG algorithm to run one sweep 
with the following modifications: 

i) For each even bond apply the local time evolution U at the bond formed by 
the free sites to 1 1//) . This is a very fast operation compared to determining the 
ground state, which is usually done instead in the finite-system algorithm. 

ii) As always, perform a DMRG truncation at each step of the finite-system 
algorithm, hence 0{L) times. 

(iii) Use White's prediction method to shift the free sites by one. 

2. In the reverse direction, apply step (i) to all odd bonds. 

3. As in standard finite-system DMRG evaluate operators when desired at the end of 
some time steps. Note that there is no need to generate these operators at all those 
time steps where no operator evaluation is desired, which will, due to the small 
Trotter time step, be the overwhelming majority of steps. 

Note that one can also perform every bond evolution operator at each half-sweep, in 
order. This does not worsen the Trotter error, since in the reverse sweep the operators 
are applied in reverse order. 

The calculation time of adaptive time-dependent DMRG scales linearly in L, as op- 
posed to the static time-dependent DMRG which does not depend on L. The diagonaliza- 
tion of the density matrices (Schmidt decomposition) scales as n^iigin^', the preparation 
of the local time evolution operator as n^^^^, but this may have to be done only rarely e.g. 
for discontinuous changes of interaction parameters. Carrying out the local time evolu- 
tion scales as n^-jgrn^; the basis transformation scales as n^-^^^m^. As m ^ Hjite typically, 
the algorithm is of order O(Ln^-jgm^) at each time step. 

The performance of this method has been tested in various applications in the context 
of ultracold atom physics but also for far-from-equilibrium dynamics 



1 27] and for spectral functions some of these applications will serve as examples in 
the following. 



FAR-FROM-EQUILIBRIUM DYNAMICS 

In this section, we consider the dynamics of a system far from equilbrium using adap- 
tive time-dependent DMRGi27l]. The following example, for which an exact solution 
is available, shows that time-dependent DMRG can also perform in situations where 
dynamical DMRG must surely fail. 

The initial state |ini) = ||...ti---i)on the one-dimensional spin- 1/2 chains is 
subjected to the dynamics of the Heisenberg model 

H = "L^lSl^i-^^K+i +JzS'nSl^, = l^K. (44) 

n n 

We set h = I, defining time to be 1/energy with the energy unit chosen as the J^y 
interaction. 

Often it is useful to map the Heisenberg model onto a model of interacting spinless 
fermions with nearest-neighbour hopping: 



H 



2' 



(45) 



In particular, the case 7- = describes free fermions on a lattice, and can be solved 
exactly. In the following we will focus on this case. Note that in that case the initial 
state with two large ferromagnetic domains separated by a domain wall in the center is 
a highly excited state; the ground state exhibits power-law decaying antiferromagnetic 
correlations. 

The time evolution delocalizes the domain wall over the entire chain; the magnetiza- 
tion profile for the initial state |ini) reads Ii28ll : 



S,M = {w{t)K\Wit)) = -1/2 £ Jjit), (46) 



where Jj is the Bessel function of the first kind, n = . . . , — 3, —2,-1,0, 1,2, 3, . . . labels 
chain sites with the convention that the first site in the right half of the chain has label 
R = 1 . As the total energy of the system is conserved, the state cannot relax to the ground 
state. The exact solution reveals a nontrivial behaviour with a complicated substructure 
in the magnetization profile, which is a good benchmark for DMRG. 

Possible errors. Two main sources of error occur in the adaptive t-DMRG: 
(i) The Trotter error due to the Trotter decomposition. For an nth-order Trotter decom- 



position ll25h . the error made in one time step dt is of order Ldt'^^^ . To reach a given 



time t one has to perform t /dt time-steps, such that in the worst case the error grows 



linearly in time t and the resulting error is of order L{dt)"t. 

(ii) The DMRG truncation error due to the representation of the time-evolving quantum 
state in reduced (albeit "optimally" chosen) Hilbert spaces and to the repeated trans- 
formations between different truncated basis sets. While the truncation error £ that sets 
the scale of the error of the wave function and operators is typically very small, here it 
will strongly accumulate as 0{Lt/dt) truncations are carried out up to time t. This is 
because the truncated DMRG wave function has norm less than one and is renormal- 
ized at each truncation by a factor of (1 — e)^^ > 1. Truncation errors should therefore 
accumulate roughly exponentially with an exponent of eLt/dt, such that eventually the 
adaptive t-DMRG will break down at too long times. The accumulated truncation error 
should decrease considerably with an increasing number of kept DMRG states m. For a 
fixed time t, it should decrease as the Trotter time step dt is increased, as the number of 
truncations decreases with the number of time steps t /dt. 

At this point, it is worthwhile to mention that our subsequent error analysis should also 
be pertinent to the very closely related time-evolution algorithm introduced by Verstraete 
et q/. iEqIi . which also involves both Trotter and truncation errors. 

We remind the reader that no error is encountered in the application of the local time 
evolution operator Un to the state 1 1//) . 

Error analysis. We use two main measures for the error: 
(i) As a measure for the overall error we consider the magnetization deviation, the 
maximum deviation of the local magnetization found by DMRG from the exact result, 

err(0 =max„|(5^Dj^G(f)) - (5^,exact(0) I- (47) 



(ii) As a measure which excludes the Trotter error we use the forth-back deviation 
FB{t), which we define as the deviation between the initial state |ini) and the state 
\fb{t)) = U {—t)U {t)\mi) , i.e. the state obtained by evolving |ini) to some time t and 
then back to t = again. If we Trotter-decompose the time evolution operator U {—t) 
into odd and even bonds in the reverse order of the decomposition of U (t), the identity 
U{-t) = U{t)'^ holds without any Trotter error, and the forth-back deviation has the 
appealing property to capture the truncation error only. 

As the DMRG setup used in this particular calculation did not allow easy access to 



the fidelity \{im\ fb{t))\ (a calculation which is not a problem in principle, see a3'2\ ). 
the forth-back deviation was defined to be the L2 measure for the difference of the 
magnetization profiles of |ini) and \ fb{t)), 

m) = (£((ini|5^|ini) - {fb{t)\Sl\fb{t))f^ (48) 

In order to control Trotter and truncation error, two DMRG control parameters are 
available, the number of DMRG states m and the Trotter time step dt. 

The dependence on dt is twofold: on the one hand, decreasing dt reduces the Trotter 
error by some power of dt" exactly as in QMC; on the other hand, the number of 
truncations increases, such that the truncation error is enhanced. It is therefore not a 
good strategy to choose dt as small as possible. The truncation error can however be 
decreased by increasing m. 
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FIGURE 8. Magnetization deviation err(f ) as a function of time for different numbers m of DMRG 
states. The Trotter time interval is fixed at dt — 0.05. Again, two regimes can be distinguished: For 
early times, for which the Trotter error dominates, the error is slowly growing (essentially linearly) and 
independent of m (regime A); for later times, the error is entirely given by the truncation error, which 
is m-dependent and growing fast (almost exponential up to some saturation; regime B). The transition 
between the two regimes occurs at a well-defined "runaway time" Ir (small squares). The inset shows a 
monotonic, roughly linear dependence of f^; on m. From Ii27il . 



Consider the dependence of the magnetization deviation err(?) on the number m of 
DMRG states. In Fig. |H1 err{t) is plotted for a fixed Trotter time step dt = 0.05 and 
different values of m. One sees that a m-dependent "runaway time" tR separates two 
regimes: for t < (regime A), the deviation grows essentially linearly in time and is 
independent of m, for t>tR (regime B), it suddenly starts to grow more rapidly than any 
power-law as expected of the truncation error. In the inset of Fig.[8l is seen to increase 
roughly linearly with growing m. As m ^ oo corresponds to the complete absence of the 
truncation error, the m-independent bottom curve of Fig. [His a measure for the deviation 
due to the Trotter error alone and the runaway time can be read off very precisely as the 
moment in time when the truncation error starts to dominate. 

That the crossover from a dominating Trotter error at short times and a dominating 
truncation error at long times is so sharp may seem surprising at first, but can be 
explained easily by observing that the Trotter error grows only linearly in time, but 
the accumulated truncation error grows almost exponentially in time. 

To see that nothing special is happening at tR, consider also Fig. 13 where the Trotter- 
error free FB{t) is plotted as a function of m, for t = 30 and t = 50. An approximately 
exponential increase of the accuracy of the method with growing m is observed for a 
fixed time. Our numerical results that indicate a roughly linear time-dependence of tR 
on m (inset of Fig. [U) are the consequence of some balancing of very fast growth of 
precision with m and decay of precision with t. 
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FIGURE 9. The forth-back error FB{t) for f = 50 and f = 30, as function of m. Here, L = 100, 
i/f = 0.05. From 1^. 

The runaway time thus indicates an imminent breakdown of the method and is a good, 
albeit very conservative measure of available simulation times. We expect the above er- 
ror analysis for the adaptive t-DMRG to be generic for other models. The truncation 
error will remain also in approaches that dispose of the Trotter error; maximally reach- 
able simulation times should therefore be roughly the similar. Even if for high precision 
calculation the Trotter error may dominate for a long time, in the long run it is always 
the truncation error that causes the breakdown of the method at some point in time. 



FINITE TEMPERATURE 

After the previous discussion on the difficulties of simulating the time-evolution of pure 
states in subsets of large Hilbert spaces it may seem that the time-evolution of mixed 
states (density matrices) is completely out of reach. It is however easy to see that a 
thermal density matrix pp = exp[— /3^] can be constructed as a pure state in an enlarged 
Hilbert space and that Hamiltonian dynamics of the density matrix can be calculated 
considering just this pure state (dissipative dynamics being more complicated). In the 
DMRG context, this has first been pointed out by Verstraete, Garcia- RipoU and Cirac[2^ 
and Zwolak and Vidal[33], using essentially information-theoretical language; it has also 
been used previously in pure statistical physics language in e.g. high-temperature series 
expansions fsT]. 

To this end, consider the completely mixed state pQ = \. Let us assume that the 
dimension of the local physical state space {\(7i)} of a physical site is n. Introduce now 



a local auxiliary state space { | T,) } of the same dimension n on an auxiliary site. The 
local physical site is thus replaced by a rung of two sites, and a one-dimensional chain 
by a two-leg ladder of physical and auxiliary sites on top and bottom rungs. Prepare now 
each rung i in the Bell state 



C7/ = T; 



(49) 



Other choices of | V^q) are equally feasible, as long as they maintain in their product states 
maximal entanglement between physical states |a,) and auxiliary states |t,;). Evaluating 
now the expectation value of some local operator O'^y acting on the physical state space 
with respect to | i/Zq), one finds 

{¥'o\0'a\¥'o) = I I -[{atri\d',®V,\GlTl)]. 
The double sum collapses to 

(V^^|0'JV/^) = ^I;(CT,-|0',|CT,), 

CT, = 1 

and we see that the expectation value of O^^ with respect to the pure state | i/Zq) living on 
the product of physical and auxiliary space is identical to the expectation value of O^^ 
with respect to the completely mixed local physical state, or 

{&^)=TXaPW^ (50) 

where 

p^ = Tr,|v/^)(^|. (51) 
This generalizes from rung to ladder using the density operator 

po = Tr^|V/o)(VA)|, (52) 

where 

L 

IV^o)=niV^o) (53) 

is the product of all local Bell states, and the conversion from ficticious pure state to 
physical mixed state is achieved by tracing out all auxiliary degrees of freedom. 
At finite temperatures /3 > one uses 

P/3 = 1 = Tr,.-/^^/2|v/o)(V^ok-^^/^ 

where we have used Eq. (l52t and the observation that the trace can be pulled out as it 
acts on the auxiliary space and e~^^/^ on the physical space. Hence, 



P/3 =Tr,|v^^)(V^^ 



(54) 



where \ = e ^^/^\\j/o). Similarly, this finite-temperature density matrix can now be 

evolved in time by considering \ Mfp{t)) = e^'^'l V^j3(0)) and pp{t) = TrT|v^/3(^))(V^j3(^)l• 
The calculation of the finite-temperature time-dependent properties of, say, a Hubbard 
chain, therefore corresponds to the imaginary-time and real-time evolution of a Hubbard 
ladder prepared to be in a product of special rung states. Time evolutions generated 
by Hamiltonians act on the physical leg of the ladder only. As for the evaluation of 
expectation values both local and auxiliary degrees of freedom are traced on the same 
footing, the distinction can be completely dropped but for the time-evolution itself. 
Code-reusage is thus almost trivial. Note also that the initial infinite-temperature pure 
state needs only m— I block states to be described exactly in DMRG as it is a product 
state of single local states. Imaginary-time evolution (lowering the temperature) will 
introduce entanglement such that to maintain some desired DMRG precision m will 
have to be increased. 



TIME-STEP TARGETTING 

The Trotter based methods for time evolution discussed above, while very fast, have two 
notable weaknesses: first, there is an error proportional to the time step T squared. This 
error is usually tolerable and can be reduced to neglible levels by using higher order 
Trotter decompositions|9]. More importantly, they are limited to systems with nearest 
neighbor interactions on a single chain. This limitation is more difficult to deal with. In 
the case of narrow ladders with nearest-neighbor interactions, one can avoid the problem 
by lumping all sites in a rung into a single supersite. Another approach would be to 
use a superblock configuration with, say, three center sites, which would allow one to 
treat two-leg ladders without using supersites. Unfortunately, these approaches become 
very inefficient for wider ladders, and are not applicable at all to general long-range 
interaction terms. 

The time-step targetted (TST) method does not have these limitations. The main 
idea is to produce a basis which targets the states needed to represent one small but 
finite time step. Once this basis is complete enough, the time step is taken and the 
algorithm proceeds to the next time step. This targetting is intermediate to previous 
approaches: the Trotter methods target precisely one instant in time at any DMRG step, 
while Luo, Xiang, and Wang's approach ll20ll targetted the entire range of time to be 
studied. Targetting a wider range of time requires more density matrix eigenstates be 
kept, slowing the calculation. By targetting only a small interval of time, a smaller price 
is paid relative to the most efficient Trotter methods. In exchange for the modest loss 
of efficiency, we gain the ability to treat longer range interactions, ladder systems, and 
narrow two-dimensional strips. In addition, the error from a finite time step is greatly 
reduced relative to the second order Trotter method. 

The procedure of Luo, et. al. for targetting an interval of time is nearly ideal: one 
divides the interval into n small steps of length e, and targets = 0), = e), = 
2e), . . ., = ne), simultaneously. By targetting these wavefunctions simultaneously, 
any linear combination of them is also included in the basis. This means than the 
basis is able describe an n -I- 1-th order interpolation through these points, making it 



for reasonable e and n essentially complete over the time interval. In the TST method 
the interval is short and n is fairly small: in the implementation of |9], n = 3 and the time 
step is similar in size to the Trotter step T, say ~ 7/10 for a spin chain. 

The Runge-Kutta (R-K) implementation of this approach is defined as follows: one 
takes a tentative time step at each DMRG step, the purpose of which is to generate a 
good basis. The standard fourth order R-K algorithm is used. This is defined in terms of 
a set of four vectors: 

l^i) = TH{t)\xif{t)), 

\k2) = TH{t + T/2)[\iir{t)) + l/2\h)], 

\h) = THit + T/2)[Mt)) + l/2\k2)]. 

\k4) = TH{t + t)[\x{r{t)) + \k3)], (55) 
where H (t) = H (t) — Eq. The state at time t + T is given by 

\yf(t + T))^^ [\ki)+2\k2) +2\k3) + \k4)] + 0{x^). (56) 

We target the state at times t, t + t/3, t + 2x/3 and t + t. The R-K vectors have been 
chosen to minimize the error in \ \(f{t + t)), but they can also be used to generate at 
other times. The states at times t + T/3 and ? + 2t/3 can be approximated, with an error 
0(T4),as 

\y/{t + T/3)) ^ Mt)) + 

+ y^[31|^i) + 141^2) + 141-^3) -5|;t4)], 
|V/(r + 2T/3)) ^ Mt)) + 

+ ^[16|fci) + 201^2) + 201^3) -21-^4)]. (57) 

Each half-sweep corresponds to one time step. At each step of the half-sweep, one 
calculates the R-K vectors (l55l) . but without advancing in time. The density matrix is 
then obtained with the target states |v^(0), IV^(^ + t/3)), \\ir{t + 2t;/ 3)), and |i/a(? + t)). 
Advancing in time is done on the last step of a half-sweep. However, we may choose to 
advance in time only every other half-sweep, or only after several half-sweeps, in order 
to make sure the basis adequately represents the time-step. For the systems of Ref. ||^, 
one half-sweep was adequate and the most efficient. The method used to advance in time 
in the last step need not be the R-K method used in the previous tentative steps. In fact, 
the computation time involved in the last step of a sweep is typically miniscule, so a 
more accurate procedure is warranted. A simple way to do this is to perform, say, 10 R- 
K iterations with step r/ 10. The relative weights of the states targetted can be optimized. 
An equal weighting is not optimal; the initial time and final time are more important. In 
Ref. [9], it was found that giving a weight of 1/3 for the first and final states, and 1/6 
for the two intermediate states, gave excellent results. 

Both the Trotter method and the TST method give very accurate results. In Fig. [121 
and Fig. [TTl we show a comparison of the methods. On a large scale, we cannot 
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FIGURE 10. The value of (5|-(16,f)5+(16,0)) computed for a 31 site S = 1 Heisenberg chain, 
computed three different times. Here the curves labeled Runge-Kutta are the TST method, implemented 
using Runge Kutta. The time step is T. The difference in results are not visible on this scale. 
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FIGURE 11. Same as for Fig.no] but showing only a small region so the differences become apparent. 



see any difference between the methods for times out to f ~ 10. If we zoom in on 
a particular region, we see the effects of the finite Trotter decomposition error, here 
falling as T^. We kept m = 300 states for the TST method, and m = 200 states for the 
Trotter methods. Typically, one finds that more states must be targetted for the TST 
method, because the targetting is over a finite interval of time rather than one instant. The 
Trotter decomposition error can be eliminated almost completely by using a higher order 
decomposition. In this case, the smaller value of m still works as well as in the lower 
order methods. This combination gives the best combination of speed and accuracy. 



SPECTRAL FUNCTIONS 

Using either the Trotter or TST methods, it is straightforward to obtain spectral func- 
tions. Typically, we are interested in the Fourier transform of a time dependent correla- 





tion function 




(EDI)) as 



(58) 



C{t) = {^\Bexp[-it{H -Eg)]A\^) 



(59) 



where Eq is the ground state energy. For evaluating this expression, we proceed as 
follows: we first obtain the ground state using the standard DMRG algorithm. We 
then apply the operator A to obtain = 0)), and evolve \^^{t)) in time, using the 
Hamiltonian with the ground state energy subtracted off. During this time evolution, we 
target both \^f{t)) and We then obtain C{t), at each time step, as 



By targetting both and we ensure that this matrix element can be obtained 

with an accuracy controlled by the truncation error. 

In forming = 0)), we use a complete half-sweep to apply A to In particular, 
if A is a sum of terms Aj over a number of sites, then we apply an Aj only when j is one 
of the two central, untruncated sites. Thus the basis is automatically suitable for Aj|0). 

During this buildup of A at step j we target both 1 0) and 1 0) • At the end of the 

sweep, we turn on the time evolution. 

For translationally invariant systems it is particularly convenient to let A and B be 
on-site operators, for example A = 5'+(j), where j is in the center of the chain. Since the 
time evolution does not evolve B, a whole set of 5's can be utilized, one for each site of 
the system, for example B = S (i) . One measurement of G{£ — j, t) can be made on each 
step of each sweep, where i is one of the two center sites with untruncated bases, so that 
no extra operator matrices need be kept to reproduce B. In this way, one simulation yields 
G{i — J, t) for a wide range of values ofi — j and t. By Fourier transforming in both space 
and time, one obtains the full spectral function for all frequency and momenta, in one 
simulation. This is in stark contrast to the most accurate frequency methods, in which 
one k and a small range of (0 are obtained in one run. 

As an example we return to the isotropic spin-1 Heisenberg chain, with the exchange 
coupling J set to unity, and A = S^{j), as above. Note that the application of 5+(j) 
constructs a localized wavepacket consisting of all wavevectors. This packet spreads out 
as time progresses, with different components moving at different speeds. The speed of 
a component is its group velocity, determined as the slope of the dispersion curve at k. 
In FigUJlwe show the local magnetization (v^(?) !//(?)) for a chain of length L = 200, 
with timestep T = 0. 1 . At ? = 0, the wavepacket has a finite extent, with size given by 
the spin-spin correlation length ^ . At later times, the different speeds of the different 
components give the irregular oscillations in the center of the packet. We kept m = 150 
states per block, giving a truncation error of about 6 x 10^^. When the wavefront reaches 
the edge of the system, we stop the simulation. Our results up to that point in time have 
very minimal finite size effects, dying off exponentially from the edges. The correlation 
function is exponentially small for \i — j\ greater than vt, with v the maximum group 
velocity. Because of this, we can specify the momentum precisely and arbitrarily, i.e. as 



c{t) = {m¥it)). 



(60) 
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FIGURE 12. Time evolution of the local magnetization {S' {x)) of a 200 site spin-1 Heisenberg chain 
after 5+ (100) is applied. From [2J. 



if the system were infinite. The broadening from having a finite system appears only in 
frequency, not momentum. 

The spectral function is defined as —1/ ;rImG(fc, ft)), where 

G{x,t) = -iC{x,t) = -/((^)|r[57(0S+(0)]|<^)) (61) 

Since G(x, t) is even in x and t, the Fourier transform is 

/"OO 

G{k,co) = 2 dt cos cotYcoskxGlxj) (62) 

In this expression, it is the real part of C{xj) that determines the spectral function; the 
imaginary part is thrown away. For an excited state with energy A above the ground 
state, this spectral function gives peaks at ±A. Alternatively, we can define the spectral 
function to have only the +A peaks. This spectral function comes from a Fourier 
transform utilizing both the real and imaginary parts of C{x,t), and both positive and 
negative times: 

1 f°° 

A{k,(o) = — dte''^'Y,^oskx{^\S-{t)S^\^). (63) 

J— oo ^ 

In this case, we utilize C{x,t) = C(a:, —t)* to obtain the negative time data. We prefer 
this latter spectral function, since it utilizes the imaginary data, but the differences are 
not large and we have not studied them carefully. 

We approximate the time integral utilizing a windowing function W{t) which goes to 
zero as ? — i> r, 

W{t). (64) 



A set of windowing functions with a number of nice properties is 
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W„(0=cos(— )". (65) 
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FIGURE 13. The single magnon spectrum of the spin-1 Heisenberg antiferromagnetic chain, for a 
system of L = 600 sites, using the Trotter method with a time step of T = 0.1, running for T ~ 100, and 
keeping up to m — 600 states, at momentum k=n. The main peak has a height of 83 at o = 0.415, close 
to the true Haldane gap of A = 0.41050(2). 1 35 1 The sharp oscillations around it are the result of numerical 
errors and windowing. The three-magnon continuum is visible, beginning at 3A. js^l 
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FIGURE 14. Same as for Fig.[T3] but at ;t = k/IQ 



These functions approach Gaussians as n ^ oo, but the function and n—\ derivatives 
vanish at ? = ±.T . If one sets Wn{t) to zero for \t\ > T, and Fourier transforms, one 
obtains a nearly Gaussian lineshape with oscillating tails falling off as co^^"^^\ We 
have used n = 4, for which the lineshape in CO has negative regions in the tails of very 
small amplitude, less than half a percent of the peak height. Another good choice is 
n = 3, giving a somewhat narrower peak at the expense of more negativity. Note that if 
the true spectral function has an isolated delta function peak, the windowed spectrum 
will have a broadened peak centered precisely at the same frequency. Thus it is possible 
to locate the single magnon line with an accuracy much better than l/T. If a continuum 
is also present nearby, the peak is less well determined. In the case of the 5=1 chain, 
for k near n the peak is isolated, but at some point near k = 0.25n the peak enters the 
two magnon continuum and develops a finite width. 



In FiglOlwe show the resuking spectrum for k= n. The resuks for the three magnon 
continuum are impressive, as the single magnon line is much larger in amplitude. In 
Figdlwe show results for k = n/lO. For this momentum, the single magnon line lies 
within the two magnon continuum, altering its shape dramatically. This shape has been 
calculated approximately using the nonlinear sigma model [13611 : the results shown are in 
good qualitatively agreement with these analytic results. 



CONCLUSION 

While the invention of efficient time-dependent DMRG methods is at the time of writ- 
ing only one and a half years old, the results achieved so far have already been impres- 
sive, indicating that the problem of highly precise time-evolutions for one-dimensional 
strongly correlated quantum is for the first time under very good control. The available 
variants cover a wide range of physical problems, with Trotter-based methods most ef- 
ficient for short-ranged Hamiltonians, and with time-step targetting methods superior 
for longer-ranged interactions. They also provide powerful alternatives for the calcula- 
tion of spectral functions in the momentum-frequency range. A nice feature is provided 
by their easy implementation in the framework of existing static finite-system DMRG 
codes. As control of quantum systems is improving experimentally, we expect the range 
of physical applications to grow strongly in the very near future. 
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